*-------------------------------------------------------------------------------
*Project: School Closures
*This code: Submission Package - Replication
*Author: Anelli, Galasso, Griselda
*-------------------------------------------------------------------------------
*Date created 14/05/22
*-------------------------------------------------------------------------------

clear all
cd "" // insert here the right working directory

* Load dataset

use "school_closure_replication.dta", clear

*************************
**# CLEAN RESPONSE DATA *
*************************

* 1. Drop if the response to all activities is 0

egen tot_act=rsum(act_class act_hmwk_superv act_hmwk_alone act_extracurric act_other) if inrange(schoolkid,2,4)
drop if act_class==0 & act_hmwk_alone==0 & act_hmwk_superv==0 & act_extracurric==0 & act_other==0

* 2. Drop if they allocate fewer than 5 hours in total

drop if tot_act<5

* 3. Drop if they allocate 0 hours to both in class and online learning

gen act_class_bef=act_class if pandemic==0
gen act_class_dur=act_class if pandemic==1

bysort ID12 (pandemic) : replace act_class_bef = act_class_bef[_n-1] if missing(act_class_bef) 
bysort ID12 (pandemic) : replace act_class_dur = act_class_dur[_n+1] if missing(act_class_dur) 

gen 	double_zero=. 
replace double_zero=1 if act_class_bef==0 & act_class_dur==0
replace double_zero=0 if inrange(act_class_bef,1,10) | inrange(act_class_dur,1,10)
drop 	if double_zero==1

* 4. Drop if they allocate 10 hours to both in class and online learning

gen 	double_ten=.
replace double_ten=1 if act_class_bef==10 & act_class_dur==10
replace double_ten=0 if inrange(act_class_bef,0,9) | inrange(act_class_dur,0,9)
drop 	if double_ten==1

* 5. Drop if they allocate 10 hours to remote learning

drop if act_class==10 & pandemic==1
	
drop act_class_bef act_class_dur


*************************************************
** GENERATE USEFUL VARIABLES - DEFINE CONTROLS *
*************************************************
gen work_outside= (job_working_curr ==1) if job_working_curr !=.
	la define job_working_curr_fromhome_d_l 0 "Parents do not work from home" 1 "Parents work from home"
	la values job_working_curr_fromhome_d job_working_curr_fromhome_d_l
	
* Create variable of time set
sort ID12 pandemic
encode ID12, gen(ID12_en)
xtset ID12_en pandemic
gen homework=act_hmwk_superv+act_hmwk_alone

gen delta_SH=d.act_hmwk_superv
gen delta_UH=d.act_hmwk_alone
gen delta_homework=d.homework
gen delta_Extra=d.act_extracurric
gen delta_Other=d.act_other

* Manually add data on returns to education
gen 	beta_SL=	0.00232 // Lavy 2015
gen 	beta_RL=	beta_SL*0.943 // Alpert et al 2016
gen 	beta_SH=0.00279 // Resnjanskij 2021
gen 	beta_UH=	0.00171 // Eren Anderson 2011

************
** DESCRIPTIVE GRAPHS *
************
do replication_graphs.do 

***********
** TABLE *
***********

do replication_table.do

***********
** MODEL *
***********

/* Figure_5: Performance Loss Due to School Closure in 2020 */

collapse (mean) mean_act_class= act_class_ mean_act_hmwk_superv= act_hmwk_superv  mean_act_hmwk_alone= act_hmwk_alone mean_act_extracurric= act_extracurric mean_act_other=act_other beta_SL beta_SH beta_UH lost_day_2020 lost_day_2021, by(cntr pandemic schoolkid)  

sort cntr pandemic
order cntr pandemic

gen D_SL = -mean_act_class if pandemic==0 
gen D_RL =  mean_act_class if pandemic==1 

bysort cntr schoolkid (pandemic): gen D_SH =  mean_act_hmwk_superv - mean_act_hmwk_superv[_n-1]

bysort cntr schoolkid (pandemic): gen D_UH =  mean_act_hmwk_alone - mean_act_hmwk_alone[_n-1]

* Sweden is a particular case where primary and lower secondary schools never
*closed.
bysort cntr schoolkid (pandemic): gen D_SL_sw =  mean_act_class - mean_act_class[_n-1]
replace D_SL=D_SL_sw 	if cntr==8 & inrange(schoolkid,2,3)  
replace D_RL=0 			if cntr==8 & inrange(schoolkid,2,3)

collapse mean_act_class mean_act_hmwk_superv mean_act_hmwk_alon mean_act_extracurric mean_act_other beta_SL beta_SH beta_UH lost_day_2020 lost_day_2021 D_SL D_RL D_SH D_UH, by(cntr pandemic)  

foreach var in  D_SL D_RL D_SH D_UH {
	bysort cntr (pandemic) : replace `var' = `var'[_n-1] if missing(`var') 
	bysort cntr (pandemic) : replace `var' = `var'[_n+1] if missing(`var') 
}

collapse  D_SL D_RL D_SH D_UH beta_SL beta_SH beta_UH lost_day_2020 lost_day_2021, by(cntr)
order cntr D_SL D_RL D_SH D_UH

* Make day of school closure consistent
replace lost_day_2020=61 if cntr==1
replace lost_day_2020=37 if cntr==3
replace lost_day_2020=24 if cntr==4
replace lost_day_2020=78 if cntr==5
replace lost_day_2020=24 if cntr==6
replace lost_day_2020=93 if cntr==7
replace lost_day_2020=19 if cntr==8
replace lost_day_2020=39 if cntr==9

* Varying Remote learning return
gen 	beta_RL_1=	beta_SL*0.943 // Alpert et al 2016
gen 	beta_RL_6= 	beta_SL*0.799144779 // rate such that total loss in Italy=-0.19SD as in the literature


foreach i in 1 6{
	gen school_loss_daily_`i'= beta_SL*D_SL + beta_RL_`i' * D_RL
	la variable school_loss_daily_`i' "Daily School Loss in Cognitive Abilities (SD)"	
	gen school_loss_2020_`i' =school_loss_daily_`i' * lost_day_2020

}

gen parental_comp_daily=beta_SH*D_SH + beta_UH* D_UH
	la variable parental_comp_daily "Parental and Student Compensation in Cognitive Abilities (SD)"
gen parental_comp_2020	=parental_comp_daily* lost_day_2020
	la variable parental_comp_2020  "Parental and Student Compensation in Cognitive Abilities (SD) in 2020"

foreach i in 1 6{
	gen delta_HC_daily_`i'=school_loss_daily_`i'+parental_comp_daily
	la variable delta_HC_daily_`i' "Overall Daily Loss in Cognitive Abilities (SD)"	
	gen delta_HC_2020_`i'		=delta_HC_daily_`i'		* lost_day_2020

}	

// compute the upper and lower bound for the graph
gen school_loss_2020_lb=school_loss_2020_1
gen school_loss_2020_ub=school_loss_2020_6
gen school_loss_2020=school_loss_2020_1
	la variable school_loss_2020 	"School Loss in Cognitive Abilities (SD) in 2020"

gen delta_HC_2020_lb=delta_HC_2020_1
gen delta_HC_2020_ub=delta_HC_2020_6
gen delta_HC_2020=delta_HC_2020_1
		la variable delta_HC_2020		"Overall Loss in Cognitive Abilities (SD) in 2020"
		
expand 4 // for graph purpose 

gen order=.
replace order=1 if cntr==8
replace order=2 if cntr==6
replace order=3 if cntr==4
replace order=4 if cntr==3
replace order=5 if cntr==9
replace order=6 if cntr==1
replace order=7 if cntr==5
replace order=8 if cntr==7

sort order
gen n=_n

gen school_loss_2020_		=school_loss_2020 		if n==1 | n==5 | n==9 | n==13 | n==17 | n==21 | n==25 | n==29
gen school_loss_2020_lb_	=school_loss_2020_lb 	if n==1 | n==5 | n==9 | n==13 | n==17 | n==21 | n==25 | n==29
gen school_loss_2020_ub_	=school_loss_2020_ub 	if n==1 | n==5 | n==9 | n==13 | n==17 | n==21 | n==25 | n==29

gen parental_comp_2020_		=parental_comp_2020		if n==2 | n==6 | n==10 | n==14 | n==18 | n==22 | n==26 | n==30	

gen delta_HC_2020_ 			=delta_HC_2020  	if n==3 | n==7 | n==11 | n==15 | n==19 | n==23 | n==27 | n==31
gen delta_HC_2020_lb_ 		=delta_HC_2020_lb  	if n==3 | n==7 | n==11 | n==15 | n==19 | n==23 | n==27 | n==31
gen delta_HC_2020_ub_ 		=delta_HC_2020_ub  	if n==3 | n==7 | n==11 | n==15 | n==19 | n==23 | n==27 | n==31

	la define n_ll 2 "SE" 6 "NZ" 10 "GER" 14 "FRA" 18 "UK" 22 "AUT" 26 "ITA" 30 "PL"
	la values n n_ll
	
gen prop_parental_compensation=	(parental_comp_2020/-school_loss_2020_1)
gen parental_comp_2020_h=parental_comp_2020_/2

tab cntr, sum(delta_HC_2020_lb_)
tab cntr, sum(delta_HC_2020_ub_)

drop if n>33

graph twoway (bar school_loss_2020_ n,   barw(1) fcolor(navy) lcolor(navy) ) ///
			 (rbar school_loss_2020_ub_ school_loss_2020_lb_ n, sort lcolor(navy%20) fcolor(navy%30)) ///0
			 (scatter school_loss_2020_lb_ n, sort mlabel(school_loss_2020_lb_) mcolor(navy) msize(2-pt) msymbol(none) ///
			 mlabsize(tiny) mlabformat(%3.2f) mlabcolor(navy) mlabgap(tiny)) ///
			 (scatter school_loss_2020_ub_ n, sort mlabel(school_loss_2020_ub_) mcolor(navy) msize(2-pt) msymbol(none) ///
			 mlabsize(tiny) mlabformat(%3.2f) mlabcolor(navy) mlabgap(tiny)) ///
			 (bar parental_comp_2020_ n, barw(1) fcolor(orange) lcolor(orange)) ///
			 (scatter parental_comp_2020_ n, msize(zero) mlabcolor(maroon) mlabel(parental_comp_2020_) mlabposition(12) mlabsize(tiny) mlabformat(%3.2f)) ///
			 (bar delta_HC_2020_ n,      barw(1) fcolor(green) lcolor(green)) ///
			 (rbar delta_HC_2020_ub_ delta_HC_2020_lb_ n, sort lcolor(green%20) fcolor(green%30)) ///
			 (scatter delta_HC_2020_lb_ n, sort mlabel(delta_HC_2020_lb_) mcolor(green) msize(2-pt) msymbol(none) ///
			 mlabsize(tiny) mlabformat(%3.2f) mlabcolor(navy) mlabgap(tiny)) ///
			 (scatter delta_HC_2020_ub_ n, sort mlabel(delta_HC_2020_ub_) mcolor(green) msize(2-pt) msymbol(none) ///
			 mlabsize(tiny) mlabformat(%3.2f) mlabcolor(navy) mlabgap(tiny)),  ///
			 ytitle("Learning Loss (SD)") ///
			 legend(order(1 "Change in cognitive abilities (SD) in 2020 due to change in school inputs" ///
						  5 "Change in cognitive abilities (SD) in 2020 due to change in parental and students' inputs" ///
						  7 "Overall change in cognitive abilities (SD) in 2020") rows(3) size(vsmall)) ///
			xlabel(2 6 10 14 18 22 26 30, valuelabel) xtitle("") xscale(range(0 33)) ///
			graphregion(color(white))	
						
gr_edit xaxis1.style.editstyle majorstyle(tickstyle(textstyle(size(small)))) editcopy
gr_edit xaxis1.major.num_rule_ticks = 0
gr_edit xaxis1.edit_tick 9 2 `"19 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 10 6 `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 11 10 `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 12 14 `"37 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 13 18 `"39 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 14 22 `"61 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 15 26 `"78 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 16 30 `"93 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
graph export "Figure_5.pdf", as(pdf) name("Graph") replace

/* Figure_6: Hours of Supervised Homework or In-Person Teaching Required to Compensate the
Learning Losses */

tab n, sum(delta_HC_2020_)
tab cntr, sum(delta_HC_2020_)
tab cntr, sum(delta_HC_2020_lb)
tab cntr, sum(delta_HC_2020_ub)

* The following calculations compute the number of days of school that would be 
* needed to compensate the loss. For example, AUT had a loss of 0.14SD; 
* since each hour of in-person teaching increases the cognitive ability by 
* 0.00232 (as per Lavy 2015), and because each day of school has about 5 hours 
* of in-person teaching in order to compensate the loss:
gen compensation_hours_=-delta_HC_2020_/0.00232
gen compensation_hours_lb_=-delta_HC_2020_lb_/0.00232
gen compensation_hours_ub_=-delta_HC_2020_ub_/0.00232
tab cntr, sum(compensation_hours_)

gen compensation_tutor_hour_=		-delta_HC_2020_/0.00279
gen compensation_tutor_hour_lb=		-delta_HC_2020_lb_/0.00279
gen compensation_tutor_hour_ub=		-delta_HC_2020_ub_/0.00279
tab cntr, sum(compensation_tutor_hour_)

bysort cntr (compensation_tutor_hour_) : replace compensation_tutor_hour_ = compensation_tutor_hour_[_n-1] if missing(compensation_tutor_hour_) 
bysort cntr (compensation_tutor_hour_lb) : replace compensation_tutor_hour_lb = compensation_tutor_hour_lb[_n-1] if missing(compensation_tutor_hour_lb) 
bysort cntr (compensation_tutor_hour_ub) : replace compensation_tutor_hour_ub = compensation_tutor_hour_ub[_n-1] if missing(compensation_tutor_hour_ub) 

sort n 

replace compensation_tutor_hour_=. if !missing(school_loss_2020_) | !missing(compensation_hours_)
replace compensation_tutor_hour_lb=. if !missing(school_loss_2020_) | !missing(compensation_hours_)
replace compensation_tutor_hour_ub=. if !missing(school_loss_2020_) | !missing(compensation_hours_)
replace compensation_tutor_hour_=. if n==4 | n==8  | n==12 | n==16 | n==20 | n==24 | n==28 | n==32
replace compensation_tutor_hour_lb=. if n==4 | n==8  | n==12 | n==16 | n==20 | n==24 | n==28 | n==32
replace compensation_tutor_hour_ub=. if n==4 | n==8  | n==12 | n==16 | n==20 | n==24 | n==28 | n==32

	la define n_lll 3 "SE" 7 "NZ" 11 "GER" 15 "FRA" 19 "UK" 23 "AUT" 27 "ITA" 31 "PL"
	la values n n_lll

twoway (bar compensation_tutor_hour_ n, sort barw(1) fcolor(gray) lcolor(gray)) ///
		(rbar compensation_tutor_hour_ub compensation_tutor_hour_lb n, sort lcolor(gray%20)  fcolor(gray%20)) ///
		(scatter compensation_tutor_hour_ub n, sort mlabel(compensation_tutor_hour_ub) mcolor(gray) msize(1-pt) msymbol(none) ///
		mlabsize(vsmall) mlabformat(%3.0f) mlabcolor(gray) mlabgap(tiny) mlabposition(9)) ///
		(scatter compensation_tutor_hour_lb n, sort mlabel(compensation_tutor_hour_lb) mcolor(gray) msize(1-pt) msymbol(none) ///
		mlabsize(vsmall) mlabformat(%3.0f) mlabcolor(gray) mlabgap(tiny) mlabposition(9)) ///
		(bar compensation_hours_ n, sort barw(1) fcolor(purple) lcolor(purple)) ///
		(rbar compensation_hours_ub_ compensation_hours_lb_ n, sort lcolor(purple%20)  fcolor(purple%20)) ///
		(scatter compensation_hours_ub_ n, sort mlabel(compensation_hours_ub_) mcolor(purple) msize(2-pt) msymbol(none) ///
		mlabsize(vsmall) mlabformat(%3.0f) mlabcolor(purple) mlabgap(tiny)) ///
		(scatter compensation_hours_lb_ n, sort mlabel(compensation_hours_lb_) mcolor(purple) msize(2-pt) msymbol(none) ///
		mlabsize(vsmall) mlabformat(%3.0f) mlabcolor(purple) mlabgap(tiny)) ///
		, xlabel(3 7 11 15 19 23 27 31, valuelabel) xtitle("") ylabel(0(10)90) ///
		ytitle("Hours") ///
		legend(order(1 "Hours of Supervised Homework required to recover the learning loss due to COVID-19" 5 "Hours of in-person teaching required to recover the learning loss due to COVID-19") row(2) size(vsmall)) graphregion(color(white)) 

	gr_edit xaxis1.style.editstyle majorstyle(tickstyle(textstyle(size(small)))) editcopy
	gr_edit xaxis1.major.num_rule_ticks = 0
	gr_edit xaxis1.edit_tick 9  3 `"19 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 10 7 `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 11 11 `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 12 15 `"37 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 13 19 `"39 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 14 23 `"61 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 15 27 `"78 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
	gr_edit xaxis1.edit_tick 16 31 `"93 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )

graph export "Figure_6.pdf", as(pdf) name("Graph") replace

/* Figure_S5: Percentage Loss of Hourly Salary and Cumulative Loss */

* Merge wage data

merge m:1 cntr using "returns_to_skills.dta"
drop _merge

sort n 

replace percimpactonhourlywagelow=. 	if missing(parental_comp_2020_) 
replace percimpactonhourlywagehigh=. 	if missing(parental_comp_2020_) 
replace percimpactonhourlywagelow=-(percimpactonhourlywagelow)*100
replace percimpactonhourlywagehigh=-(percimpactonhourlywagehigh)*100

egen 	percimpactonhourlywage=rmean(percimpactonhourlywagelow percimpactonhourlywagehigh)

replace cumulatedimpactoverall_low=. 	if missing(delta_HC_2020_)
replace cumulatedimpactoverall_hig=.	if missing(delta_HC_2020_)
replace cumulatedimpactoverall_low=-cumulatedimpactoverall_low
replace cumulatedimpactoverall_hig=-cumulatedimpactoverall_hig

egen 	cumulatedimpactoverall=rmean(cumulatedimpactoverall_low cumulatedimpactoverall_hig)
	la variable cumulatedimpactoverall "Cumulative Loss ($) over 40 years of career"
	
graph twoway (bar percimpactonhourlywage n,  barw(1) fcolor(black%30) lcolor(black) ) ///
	(rcap percimpactonhourlywagelow percimpactonhourlywagehigh n, sort lcolor(black)) ///
	(scatter percimpactonhourlywagelow n, sort mlabel(percimpactonhourlywagelow) msymbol(none) ///
	mlabsize(tiny) mlabformat(%3.2f) mlabcolor(black) mlabgap(tiny) mlabposition(9)) ///
	(scatter percimpactonhourlywagehigh n, sort mlabel(percimpactonhourlywagehigh)  msymbol(none) ///
	 mlabsize(tiny) mlabformat(%3.2f) mlabcolor(black) mlabgap(tiny) mlabposition(9)) ///
	(bar cumulatedimpactoverall n, sort yaxis(2) barw(1) fcolor(forest_green%50) lcolor(forest_green)) ///
	(rcap cumulatedimpactoverall_low cumulatedimpactoverall_hig n, sort yaxis(2) lcolor(forest_green)) ///
	(scatter cumulatedimpactoverall_low n, sort yaxis(2) mlabel(cumulatedimpactoverall_low) msymbol(none) ///
	mlabsize(tiny) mlabformat(%13.0fc) mlabcolor(forest_green) mlabgap(tiny) mlabposition(3)) ///
	(scatter cumulatedimpactoverall_hig n, sort yaxis(2) mlabel(cumulatedimpactoverall_hig)  msymbol(none) ///
	mlabsize(tiny) mlabformat(%13.0fc) mlabcolor(forest_green) mlabgap(tiny) mlabposition(3)), ///
	ytitle("Percentage Loss of Hourly Salary") ///
	legend(order(1 "Percentage Loss of Hourly Salary due to COVID-19" ///
	5 "Cumulative Loss ($) over 40 years of career" ) rows(2) size(vsmall)) ///
	xlabel(2.5 6.5 10.5 14.5 18.5 22.5 26.5 30.5, valuelabel) xtitle("") ///
	graphregion(color(white)) ///
	note("To compute the percentage losses on hourly salary and the cumulated dollar loss over 40 years" "of career we use country-specific returns to skill and the median wages for each country.", size(vsmall))
			
graph twoway (bar percimpactonhourlywagelow n,  barw(1) fcolor(gray) lcolor(gray) ) ///
	(rbar percimpactonhourlywagelow percimpactonhourlywagehigh n, sort lcolor(black%20) fcolor(black%20)) ///
	(scatter percimpactonhourlywagelow n, sort mlabel(percimpactonhourlywagelow) msymbol(none) ///
	mlabsize(tiny) mlabformat(%3.2f) mlabcolor(black) mlabgap(tiny) mlabposition(9)) ///
	(scatter percimpactonhourlywagehigh n, sort mlabel(percimpactonhourlywagehigh)  msymbol(none) ///
	mlabsize(tiny) mlabformat(%3.2f) mlabcolor(black) mlabgap(tiny) mlabposition(9)) ///
	(bar cumulatedimpactoverall_low n, sort yaxis(2) barw(1) fcolor(forest_green) lcolor(forest_green)) ///
	(rbar cumulatedimpactoverall_low cumulatedimpactoverall_hig n, sort yaxis(2) lcolor(forest_green%30) fcolor(forest_green%30)) ///
	(scatter cumulatedimpactoverall_low n, sort yaxis(2) mlabel(cumulatedimpactoverall_low) msymbol(none) ///
	mlabsize(tiny) mlabformat(%13.0fc) mlabcolor(forest_green) mlabgap(tiny) mlabposition(3)) ///
	(scatter cumulatedimpactoverall_hig n, sort yaxis(2) mlabel(cumulatedimpactoverall_hig)  msymbol(none) ///
	mlabsize(tiny) mlabformat(%13.0fc) mlabcolor(forest_green) mlabgap(tiny) mlabposition(3)), ///
	ytitle("Percentage Loss of Hourly Salary") ///
	legend(order(1 "Percentage Loss of Hourly Salary due to COVID-19" ///
	5 "Cumulative Loss ($) over 40 years of career" ) rows(2) size(vsmall)) ///
	xlabel(2.5 6.5 10.5 14.5 18.5 22.5 26.5 30.5, valuelabel) xtitle("") ///
	graphregion(color(white)) 			

gr_edit yaxis1.reset_rule 0 3 0.5 , tickset(major) ruletype(range) 
gr_edit yaxis2.reset_rule 0 40000 10000 , tickset(major) ruletype(range) 
gr_edit yaxis2.title.text.Arrpush "Cumulative loss ($)"
gr_edit yaxis2.major.label_format = `"%13.0fc"'
gr_edit xaxis1.style.editstyle majorstyle(tickstyle(textstyle(size(small)))) editcopy
gr_edit xaxis1.major.num_rule_ticks = 0
gr_edit xaxis1.edit_tick 1 2.5 `"SE"', tickset(major)
gr_edit xaxis1.edit_tick 2 6.5 	"NZ"  , tickset(major)	 
gr_edit xaxis1.edit_tick 3 10.5 "GER" , tickset(major) 
gr_edit xaxis1.edit_tick 4 14.5 "FRA" , tickset(major) 
gr_edit xaxis1.edit_tick 5 18.5 "UK"  , tickset(major) 
gr_edit xaxis1.edit_tick 6 22.5 "AUT" , tickset(major) 
gr_edit xaxis1.edit_tick 7 26.5 "ITA" , tickset(major) 
gr_edit xaxis1.edit_tick 8 30.5 "PL"  , tickset(major)
gr_edit xaxis1.edit_tick 9  2.5  `"19 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 10 6.5  `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 11 10.5  `"24 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 12 14.5  `"37 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 13 18.5  `"39 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 14 22.5  `"61 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 15 26.5  `"78 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
gr_edit xaxis1.edit_tick 16 30.5  `"93 days"', tickset(major) editstyle(tickstyle(textstyle(size(small)))) editstyle(tickstyle(textgap(3)) )
graph export "Figure_S5.pdf", as(pdf) name("Graph") replace
